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ABSTRACT 


The effects of extended precision computing and other numerical techniques are 
evaluated for the Fourier matching method (FMM) acoustic scattering model, initially 
developed by Assistant Professor D. Benjamin Reeder, CDRAJSN (NFS), and Professor 
Timothy K. Stanton (MIT/WHOI). Theory on acoustic scattering, reverberation, 
scattering models, conformal mapping, scatterer boundary conditions, floating-point 
arithmetic, computational error, and extended precision computing is presented as a 
foundation for research development. The paper presents an assessment of the effects of 
numerical techniques on model output with the initial expectation of obtaining a more 
accurate, converged solution at higher frequencies, higher modal combinations, and 
greater eccentricities of scatterer shape. Comparisons to results from Reeder and Stanton 
(2004) demonstrate effects of executed techniques. Analysis includes an evaluation of 
the relationship between variable precision settings and computational time, gains in the 
useful frequency regime of the FMM, and numerical analysis benefits. Demonstrated 
techniques confirm that increased precision has a positive effect on model performance. 
The utility of other numerical techniques is discussed, and limitations of current 
computer systems and other shortfalls are illustrated. A feasibility assessment for Navy 
use of the FMM and recommendations for further improvements to the FMM are 
included. 
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I. INTRODUCTION 


A. SIMULATIONS OF ACOUSTIC SCATTERING 

The design of a eapable sonar system requires a thorough understanding of the 
physieal meehanisms related to the propagation of underwater sound. Two of these 
meehanisms, seattering and reverberation, eneompass the refleetion of sound from both 
oeeanie boundaries and objects suspended within the water column. Scattering from 
submerged objects typically dissipates energy away from the source-receiver path and 
obscures returned echoes through a phenomenon called reverberation, limiting the 
performance of sonar systems in certain ocean environments. These two phenomena, 
scattering and reverberation, can be only marginally simulated even in today’s most 
advanced sonar systems, and thus limit operational effectiveness. 

A school of fish can be considered a volume of scatterers that produces 
reverberation when insonified by an incident pulse from an active sonar system. If 
scattering by the school of fish could be modeled based upon an enhanced understanding 
of the physical mechanisms that produced the echo or reverberation, then this 
reverberation could be parameterized within the active sonar equation. The amount of 
reverberation in the equation could be reduced based upon an understanding of the origin 
and character of the associated scattering events. By incorporating more finely-tuned 
scattering simulations, next-generation sonar systems may take advantage of improved 
signal-to-noise ratios, which will result in the increased probability of detection of real 
targets. 

This research focuses on improving one such modeling approach, known as the 
Fourier matching method (FMM). The FMM is valid for all frequencies and scatterer 
shapes but is numerically limited in its utility and application to acoustic scattering by 
fairly regular shapes, because of several computational constraints. This research aims at 
addressing and conquering some of the computational limitations that are impeding 
implementation of the FMM into operational systems. 
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B, HISTORY OF SCATTERING MODELS AND THEIR APPLICATION 

The problem of simulating acoustic scattering phenomena has challenged 
acousticians and mathematicians for decades. Lord Rayleigh studied scattering from a 
simple sphere during the late 1800’s (Rayleigh, 1945), while Victor Anderson published 
a paper on the solution of sound scattering from a fluid sphere in 1950 (Anderson, 1950). 
However, complex shapes present the most difficult challenge to acousticians. Several 
numerical solutions have been developed for complex shapes, including the boundary 
element method (Tobacman, 1984; Francis, 1993), the T-matrix method (Waterman, 
1968; Varadan et al, 1982; Lakhtakia et al, 1984; Hackman and Todoroff, 1985), and 
the mode matching method (Yamashita, 1990). However, these numerical methods all 
have limitations in frequency range, useful boundary conditions, scatterer surface type, 
dimensions of the scatterer, or mathematical and computational efficiency (Reeder and 
Stanton, 2004). 

Reeder and Stanton (2004) introduced a general scattering formulation for 
calculating the far-field scattered sound pressure from irregular, axisymmetric, finite- 
length bodies. Their work incorporated a two-dimensional conformal mapping approach, 
which they adapted to scattering by finite-length bodies for three boundary conditions— 
soft, rigid, and fluid. Although their method was numerically efficient in its formulation, 
the range of parameters for which the solution converged was believed to be limited by 
certain factors, such as computer precision, which is an inherent limitation of many 
numerical models. This posed restrictions to the useful frequency range and to 
eccentricity of scatterer shape (Reeder and Stanton, 2004). 

C. IMPORTANCE - THE REVIVAL OF OCEANOGRAPHY 

The driving influence behind the choice of this research topic is the renewed 
interest in the ocean acoustics field within the United States Navy’s Naval Meteorology 
and Oceanography Command. It is imperative that the U.S. Navy improve upon its 
existing sonar systems to meet future warfighting requirements, in order to ensure 
undersea dominance. Future conflicts between the U.S. and adversaries with diesel 
submarines operating in the challenging littoral environment could demand considerable 
exploitation of underwater capabilities to maintain the acoustic advantage. Advancement 
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of critical technologies is a present day priority for the U.S. Department of Defense. 
Further development of the FMM may offer potential improvements in the eapabilities of 
active sonar systems, allowing employment of a superior anti-submarine warfare tool in 
the oeeanie battlespaee. 

Submarines eontinue to get quieter and quieter, and sonar systems must be 
improved to maintain undersea dominanee against these threats, particularly in the 
shallow water environment. Incorporation of scattering models like the FMM into sonar 
systems may someday enable threat deteetion in an intensely reverberant environment, 
whieh is very difficult with today’s systems. 

The operational environment requires timely proeessing of environmental data 
and aeoustic returns during the prosecution of an enemy threat. Time is an important 
limitation in the real world, espeeially in combat situations. If the FMM eannot be used 
in a timely fashion by operational eomputer systems, then it is still only a researeh tool. 
With this in eonsideration, an analysis of eomputational time versus preeision is 
ineorporated into the results of this paper. 

While this is still 6.1 (basic) research, the FMM will hopefully one day be 
ineorporated into an advaneed anti-submarine warfare (ASW) or undersea warfare 
(USW) sonar system. Other potential applieations for the FMM inelude mine hunting, 
port seeurity systems, monitoring of fish populations, or oeeanographie researeh. The 
FMM is a progression of seienee. With eontinued development, this teehnology will 
enable oeeanographers and operators alike to simulate the undersea environment more 
accurately. Tools sueh as the FMM eould one day help operators maintain or assert 
undersea superiority in ways whieh are yet to be eoneeived. 

D, THESIS STATEMENT 

The two-dimensional Fourier matehing method (FMM), introdueed by DiPema 
and Stanton (1994), established a eonformal mapping approaeh, whieh maps the spaee 
variables to a new eoordinate system and matehes a eonstant radial eoordinate to the 
surfaee of the seattering body. Reeder and Stanton (2004) extended this approaeh for 
axisymmetric, finite-length bodies of revolution. The extended approaeh enabled the 

simulation and prediction of scattering events by more eomplex, bounded seattering 
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surfaces. Reeder and Stanton (2004) obtained results using the Institute for Eleetrical and 
Eleetronies Engineers (IEEE) double preeision format in which the floating-point 
numbers were stored in a 64-bit word, aeeurately representing numbers to about 16 
deeimal digits. Results of the EMM at double preeision are well-known and ean be 
duplieated with the current model. 

The two primary limitations of the EMM are (1) that it requires the seattering 
body to be an axisymmetric body of revolution, and (2) eomputer precision (Reeder and 
Stanton, 2004). The limitations of the EMM bring several researeh questions to mind: 

1. In what ways ean the Reeder and Stanton (2004) results be improved? 

2. Will inereased preeision eomputing lead to an improved eonverged 
solution? 

3. Are there other potential improvements that eould help the EMM produee 
a more aeeurate deseription of seattering phenomena? 

4. What are the limiting faetors for eomputational aeeuracy and effieieney 
and what numerieal teehniques ean be applied to the EMM to eireumvent 
these limiting faetors? 

5. How do the results of the improved model eompare with previous model 
results and real data eolleeted in laboratory situations? 

The focus of this research is improving the performanee of the EMM for 
applieation to seatterers of inereased eeeentrieity at higher frequeneies. Speeifieally, this 
researeh aims to improve EMM performanee by inereasing eomputer preeision and by 
evaluating the numerieal teehniques that are used within the model in order to aehieve 
greater aeeuraey and effieieney. Inereases in preeision are assessed to determine the 
added value. Various values of precision are utilized for investigations of eomputational 
time, beeause time is an operational eonstraint. Other possible improvements to the 
Reeder and Stanton model, such as numerical techniques, are investigated and 
ineorporated into the researeh model where applieable. 

This researeh will show whether or not eomputer preeision plays an important 
role in the ealeulated representation of sound seattering events. The extent to whieh 
extended preeision may improve model output is eurrently unknown. The exaet 
relationship between eomputational precision and increased accuracy for a series solution 
is not well-understood. This study attempts to charaeterize this relationship. 
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The FMM has future applieability to active sonar systems in minimizing volume 
reverberation by subtracting the unwanted acoustic returns of recognized scatterers from 
the total reverberation due to active sonar transmissions. Subtracting these known echoes 
from the reverberation will result in an improvement in signal-to-noise-ratio. Any gains 
in the operational utility of the FMM demonstrated here will only be amplified as 
computer precision is extended on the average computer system that will be employed by 
tomorrow’s Navy. 
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II. THEORY AND FUNDAMENTALS 


In the following paragraphs, theory and background information is presented for 
the development of the research model. Because of the breadth of topics associated with 
this research, theory on acoustic scattering, reverberation, scattering mathematics and 
models, applied computer terms and processes, computational error, numerical analysis, 
and assumptions and approximations is presented. Each of these subjects was considered 
in developing the approach for improving the EMM. General scattering formulations are 
described to show exactly where research efforts are focused. 

A. ACOUSTIC SCATTERING - CHARACTERIZING THE PHENOMENON 

Acoustic signals travel extremely long distances in the ocean environment in 
comparison with other forms of energy. Eor example, sound energy has an extremely 
small rate of attenuation in the ocean in contrast with light. This makes the acoustic 
signal a very useful tool to use in the ocean for many applications. However, objects 
encountered by the propagating acoustic pressure waves cause scattering. The range of 
propagation depends upon the attenuation rate in a continuous, fluid medium and the 
presence of inhomogeneities with different material and acoustic properties which cause 
scattering. The oceans contain many types of inhomogeneities, including schools of fish, 
bubbles, debris, suspended silt, and submarines. Each of these objects has a physical 
boundary and thus “intercepts and reradiates” some acoustic energy that interacts with 
that boundary (Urick, 1983). As the incident energy is reflected and reradiated outward 
in multiple directions, it is said to be “scattered.” 

It is important to note that the scattering problem is only one problem of many 
that are associated with the sonar equations. 

The sea is a moving medium containing inhomogeneities of various kinds, 
together with irregular boundaries, one of which is in motion. Multipath 
propagation is the rule. As a result, many of the sonar parameters 
fluctuate irregularly with time, while others change because of the 
unknown changes in the equipment and the platform on which it is 
mounted. Because of these fluctuations, a “solution” of the sonar 
equations is no more than a best-guess time average of what is to be 
expected in a basically stochastic problem. (Urick, 1983) 
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Achieving an accurate scattering simulation will reduce the amount of uncertainty in a 
complex problem. Scattering events affect the transmission loss and reverberation level 
components of the active sonar equation differently and are not very-well parameterized 
in current sonar system models. These poorly described phenomena are only a minute 
part of a largely stochastic and difficult problem. An accurate scattering simulation will 
reduce the uncertainty of the reverberation problem at its source—the individual 
scatterer. 

Scattering from simple shapes, like spheres and cylinders, can be solved exactly 
through separation of variables. However, it is extremely difficult to obtain an exact 
analytical solution to scattering events from asymmetric and complex shapes. Many 
acousticians have tried to formulate solutions with limited success. 

Some typical applications of the FMM for use in operational sonar systems might 
include scattering from air bubbles or fish swim bladders, which would have similar 
acoustical properties and boundary conditions. Fish swim at various depths in the ocean, 
but most air bubbles occur due to breaking wave processes and ship propellers at 
relatively shallow depths. Large numbers of microbubbles per unit volume have 
significant effects on near-surface, or shallow water, sound propagation (Medwin and 
Clay, 1998). 

When a bubble is insonified by an incident acoustic signal, the bubble reacts with 
a compression and rarefaction in response to the sound wave. The frequency of the 
acoustic signal and the size of the bubble dictate the response of the bubble to the 
stimulus. A large fraction of the incident energy is reflected in all directions by the 
pulsating bubble and the remaining energy is converted into heat. At a certain frequency, 
a bubble of specified diameter resonates and produces maximum extinction of the 
incident sound wave by converting it into the largest possible scattered wave and heat 
energy (Urick, 1983). 

B, REVERBERATION AND THE SONAR EQUATIONS 

The cumulative sum of all scattering events from all of the scatterers within a 
given volume is called “reverberation” (Urick, 1983). Ol’Shevskii (1967) defines 

reverberation as “a process describing the time variation of the total scattered sound field 
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observed at the point of reception following transmission of a sound signal.” 
Reverberation is generally considered the primary limitation on performance of active 
sonar systems. Reverberation that takes place within the water column, aside from its 
boundaries, is called “volume reverberation” (Urick, 1983). This occurs when scattered 
sound is returned back to its source in a monostatic arrangement. Volume reverberation 
occurs in association with a field of scatterers. 

In the ocean, the scatterers that cause volume reverberation can almost never be 
directly determined. Usually a speculation has to be made as to what types of fish, 
invertebrates, bivalves, or other scatterers are present, introducing significant uncertainty 
into the reverberation level term of the active sonar equation. If the echoes from certain 
known scatterers can be accurately modeled, then real world scatterers can also be 
modeled. Information about the underwater environment can be deduced based upon 
known reverberation characteristics. 

Early investigations of reverberation as a stochastic process were conducted by 
Yu. M. Sukharevskii (Ol’Shevskii, 1967). In 1947, he published a paper which indicated 
that the reverberation process could be described to an acceptable approximation by the 
sum of a large number of elementary scattered acoustic signals (Ol’Shevskii, 1967). If 
the ocean provided homogeneous fields of scatterers, the reverberation problem would be 
easy to solve. However, reverberation is almost always caused by a heterogeneous field 
of scatterers. This poses another problem. 

It is important to understand the scattering from individual scatterers first before 
the reverberation problem can be addressed. Then the reverberation research can be 
directly applied to real world problems. Ol’Shevskii (1967) also said, “Once, however, a 
definite hypothesis is adopted regarding the distribution of scatterers in the ocean, as well 
as their possible sizes and acoustical properties, it is possible to carry out a fairly 
complete analysis of the statistical characteristics of reverberation without analyzing in 
detail the scattering by all possible types of inhomogeneity.” If individual scattering 
events from scatterers such as fish can be mastered, then there is no need to estimate 
changing fish populations for acoustic purposes, and the reverberation problem can be 
solved with more confidence. 
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This research focuses on individual scattering events. Once the physics of the 
individual events are modeled correctly and with a certain degree of real-world utility, 
then the formulation can be applied to the largely stochastic problem of volume 
reverberation. 

C. SCATTERING MATHEMATICS AND MODELS 

1, Scattering Formulation - General Solution 

The modeling of scattering phenomena begins with the three-dimensional scalar 
wave equation; 

( 1 ) W^p = {\lc^){d^pldt^), 

where p is the acoustic pressure, is the Laplacian operator, c is the speed of sound, 
and t is time. Exact analytical solutions to the wave equation necessitate that the 
scatterer’s surface coincide with the locus of all points for which the radial coordinate is a 
constant (Reeder and Stanton, 2004). This requirement can be met for simple geometries, 
such as a sphere, an infinitely long cylinder, or a prolate spheroid. For these geometries, 
scattering phenomena can be formulated mathematically with an exact answer. 

In order to better understand acoustic scattering from a sphere, the general 
solution from Victor Anderson’s paper “Sound Scattering from a Fluid Sphere” (1950) 
will be analyzed. Anderson’s paper examines scattering from a fluid spherical scatterer, 
for which viscosity and heat conduction effects are assumed to be negligible. The sphere 
has dimensions that are similar to the wavelength of the acoustic energy with which it is 
insonified. The sphere has radius a and contains fluid 1, which has density p' and sound 
velocity c'. The sphere is centered at the origin of a polar coordinate system with radius 
denoted by r, azimuthal angle (p ranging from 0 to 2n, and polar angle z? ranging from 0 
to TT. The geometry of this sphere is shown in Figure 1. Surrounding this sphere is fluid 
2, which has different acoustical properties than the sphere that are designated by p and c 
(Anderson, 1950). 
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Figure 1. Spherical geometry for the Anderson (1950) solution to scattering from 
a fluid sphere. Azimuthal angle ^ is measured from the positive jc-axis 
in the jcj-plane and ranges from 0-2;r, The polar angle is measured 
from the positive z-axis and ranges from 0-;r, The sphere has radius, a. 

Far-field scattering is assumed, giving the aeoustie signal ample time and spaee to 
spread signifieantly from its souree. Thus, the ineident aeoustie signal is insonifying the 
seatterer in straight and parallel wave fronts. This is known as the plane-wave 
approximation (Medwin and Clay, 1998). The sphere is insonified by a plane aeoustie 
wave with pressure amplitude P and angular frequeney CO traveling parallel to the polar 
axis in the z direetion. With this ehoiee of ineident wave the dependenee on angle cp is 
removed, so that only r and must be eonsidered (i.e., baekseattering from the fluid 
sphere is identieal for all ^). Onee the sphere is insonified, the internal pressure 
beeomes p’ and a seattered spherieal wave with aeoustie pressure p emanates outward 
from the sphere. The purpose of this solution, and of the FMM, is to ealeulate the 
pressure amplitude p of this seattered wave at large ranges from the sphere. 

Along the outer edge of the sphere, where r = a, the pressure and the normal 
eomponent of partiele veloeity, u, must be eontinuous. Under these eonditions, the 
following equations should be eonsidered: 

(2) p{ci) + p^ia) = p\a ), and 
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( 3 ) 


i/^(a) + Mo^(a) = w/(a), 


where is the ineident acoustie pressure and Ur is the radial component of the particle 
velocity. The radial component of particle velocity has an incident component, Mo,r, and 
an internal component, u' , while the scattered component is denoted by Ur (Anderson, 
1950). 

The solution for acoustic pressure, p, must satisfy the three-dimensional wave 
equation given by Equation (1). The solution of the scalar wave equation is assumed to 
be time harmonic, with e~'“‘ , where CO is the angular frequency. With this assumption, 
the scalar wave equation can be transformed into the three-dimensional Helmholtz 
differential equation: 

(4) V^p{x,y,z) + k^p{x,y,z) = 0, 

In this equation, k is wave number, and k = (o! c = 27t / A, where A is the wavelength of 
the acoustic energy. For the spherical coordinate system used in this formulation, the 
Helmholtz equation is: 

(5) V^p{r,i},(p) + k^p{r,i},<p) = Q. 


In order to solve this equation, it must be expanded into: 


( 6 ) 


1 d ( 2 


dr 


dr 


-I-- 


1 


r sm 




di} 


-I-- 


\ _ ^ 

r^ sinz? dcp^ 


+ k^p = 0 , 


and then separation of variables must be carried out, with: 

(7) p{r,t},<p) = R{r)©{P)^{(p), 


where the cp terms are known to be either cos m(p or sinm^ due to sinusoidal 

d^ P 2 2 

dependence. The term —^ can be replaced with -m p , where -m is a separation 

dcp 

constant (Haberman, 1998). In addition, the time harmonic portion of the scattering 
solution of unit amplitude is p = e ~'“‘, which will be seen in the solution further on. 
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With some simplification, the separated expression becomes: 


( 8 ) 


R dr 


M 

dr 


+ k^r^ = 


-1 d 


f 


0sinz? dd 


■ 

sinz/- 

V did 


m 


sin^ id 


= ^- 


where jl is another separation constant (Haberman, 1998). The left side of this equation 
can be arranged into what is known as Sturm-Liouville form: 


( 9 ) 


d ( 2 


dr 


dr 


+ {k^r^ - fj)R = 0 (Haberman, 1998). 


This ordinary differential equation is the eigenvalue problem in r which can be solved 
with spherical Bessel functions. 

The remainder of Equation (8) above is: 


( 10 ) 

which can be simplified to: 

( 11 ) 


-1 d ( . 

smz/- 

did 


0sinz? did 


m 


sin id 




did 


(• 

f 

sinzi*- 

+ 

V dd) 

V 


jLlsmd- 


m 


.2 ^ 


sint? 


0 = 0 


where 0< d <n. This is now the eigenvalue problem in d. Assume boundary conditions 
0(0) < oo {finite) , @{7r) < oo {finite) , R{a) = 0, and |f?(0)| < oo {finite) . Let 

X = cos d to simplify the equation in d. The derivatives can be handled with the chain 
rule: 


(12) -^ = -^— = -smd— (Haberman, 1998). 

dd dd dx dx 


After dividing by sin d and recognizing that sin^ d=\- cos^ d = \-x^, the above 
equation in d becomes: 


(13) 


A 

dx 


2. d@ 


(1-xO — 

dx 


+ 

J V 


A-: 


m 


2 ^ 


0 = 0 


J 


(Haberman, 1998). 


In order to get a bounded solution at x = ±1, the simplification ji = n{n +1) is used. For 
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m = 0, where there is no dependenee on (p , the differential equation is then: 


(14) — (1“-^^)- +n(n + l)0 = O (Haberman, 1998). 

dh L dx \ 

The bounded solutions to this equation are ealled Legendre polynomials, whieh satisfy 
Rodrigues’ formula: 

(15) = (Haberman, 1998). 

2 n\dx 

The first few Legendre polynomials take the following forms: 

« = 0 : Pq{x) = 1 

(16) n = \:P^{x) = x = cosd (Haberman, 1998). 

n = 2:P^{x) = ^ (3x' -1) = ^(3eos2h»+1) 


Thus, the Legendre polynomials of order n, represented byR„(eost?), are the solution in 
the d direetion for m = 0. The assoeiated Legendre functions which are incorporated 
into the FMM formulation apply when m>0 and are similar to the Legendre 
polynomials shown here for m = 0. 

Revisiting the Sturm-Liouville Equation (9) from above, with fx = n{n +1), the 
equation becomes: 

(17) —(fy—) + (kV'-n(n + l))R = 0 

dr dr 


for n>m with fixed m. The solution at r = 0 must be bounded and R{a) = 0, where a 
is the radial coordinate at the scatterer’s surface. If kr is considered a separate variable 
and^ = kr , then there are two solution forms to Equation (17). The first solution is the 
spherical Bessel function of the first kind: 


(18) 





(Abromowitz and Stegun, 1965). 
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The second solution is the spherical Bessel function of the second kind, otherwise known 
as the spherical Neumann function; 


(19) 





(O = (-!)”"' 



(O 


(Abromowitz and Stegun, 1965). 

The solutions j„{C) the solutions to the partial differential equation for r. 


The addition of the spherical Bessel function and the spherical Neumann function yields 
the spherical Hankel function of the first kind which is of the form: 


( 20 ) K"(C} = J,.(0+‘n.(C}, 


which gives spherical harmonics of outward radiating acoustic energy. This is the 
general solution for the eigenvalue problem in the radial direction. 

Putting all of the derived components together yields the general solution for the 
fluid sphere; 


( 21 ) 


P = J^APn(COSi}) 

n=0 


Jnikr) 

n„ikr) 


-iO)t 

e 


(Anderson, 1950), 


where p is the acoustic pressure, ^4^ is the summation coefficient, F^(cosz^) is the 
Legendre polynomial, j\{kr) is the spherical Bessel function, n^{kr) is the spherical 

Neumann function, and e~“^‘ is the sinusoidal time-harmonic portion of the solution. 
This equation gives spherical harmonics of acoustic pressure. 

The solution for the internal wave, where r < a , is; 


(22) p = |;s.F.(cos . 

n-0 

The scattered wave in the region where r > a is then; 

(23) p = '^ (cos P) [y; {kr) + {kr)] , 

«=0 
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which is the portion of the solution that is most applicable to this research. However, all 
parts, including the incident acoustic wave; 

(24) =P.|;(-,r(2n + l)/>,(cos!>)y,(fc-)e'-", 

n=0 

are required to solve for the coefficients and . Here, Pq is the incident acoustic 
pressure amplitude from the source. The equation for the radial component of particle 

velocity, can be used to obtain the equations for the particle velocities 

ypc)d{kr) 

of each of the three individual waves—incident, internal and scattered. Then, equations 
for the particle velocities and the acoustic pressure equations can be solved 
simultaneously to yield the following expression for ; 

(25) H„=-Po(-0"(2n + l)/(l + /C„). 

Substituting this into Equation (21) yields: 

(26) p = -Fq^ [(-/)" {2n +1)/(I + )]X(cos [j„(kr) + {kr)]e ~'“‘, 

n—0 

which is the general solution for the total acoustie pressure external to the sphere 
(Anderson, 1950). 

The Anderson (1950) solution is directly related to the math incorporated into the 
EMM, but the EMM eontains modifications for mathematical improvements like 
geometry parameterizations and changes in boundary conditions under different physical 
circumstances. In the case of baekscattering by a fluid sphere, the EMM formulation and 
predictions are identieal to Anderson (1950) (Reeder and Stanton, 2004). 

2. Fourier Matching Method 

The Eourier Matching Method (EMM) was introduced by DiPerna and Stanton 
(1994). This method introduced a conformal mapping tactic for the prediction of far-field 
acoustic scattering for the case of an infinitely long cylinder with non-circular cross 
section. The method incorporates the conformal mapping of variables to a different 
coordinate system, in which the eonstant radial coordinate exactly matches the scatterer’s 
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surface. The eonformal mapping method generates a transformed Helmholtz equation, 
whieh is solvable even for irregularly shaped objeets (Reeder and Stanton, 2004). The 
new Helmholtz equation is solved with a given set of boundary eonditions, and the 
seattering formulation yields values for baekseattering amplitude ifbs), target strength 
(IS), and redueed target strength {RTS) for a preseribed array of frequeneies. These 
parameters are diseussed along with the FMM formulation presented here. 
a. General Solution Formulation of the FMM 
The general solution ineorporated into the FMM is similar in form to the 
Anderson (1950) general solution developed above. The formulation for the FMM 
extended to axisymmetrie, finite-length bodies by Reeder and Stanton (2004) is very 
similar to that used in the two-dimensional solution developed by DiPerna and Stanton 
(1994). The FMM is formulated with spherieal wave funetions, including spherical 
Bessel funetions, spherieal Neumann funetions, spherieal Hankel funetions, and 
assoeiated Legendre funetions. This formulation, like the Anderson (1950) solution, also 
starts with the wave equation, whieh is expressed as the Helmholtz equation in spherieal 
eoordinates given in Equation (5). The general solution to this equation is represented in 
the FMM as; 

I; I; (cosClJWe'-’+ y I; 

n=-oo m=-oo «=— 

(Reeder and Stanton, 2004). 

In this equation, j\{kr) is the spherieal Bessel funetion of order n, /^'"(eos(t?)) is the 
assoeiated Legendre funetion of degree n and order m, and h^^\kr) is the spherieal 
Hankel funetion of the first kind with order n. The total pressure external to the 
seattering body is denoted by p‘^‘{r, i}, cp). The first term in the general solution equation 
above eharaeterizes the ineident aeoustie pressure, and the seeond term corresponds to 
the seattered aeoustie pressure. The seattered field eoeffieients are denoted by B„m, and 
they ean be ealeulated by using the known ineident field eoeffieients, Anm- 

\m =+ —^^-^'^E„'”(eos(z?o)) (Reeder and Stanton, 2004). 

T{n + m + \) 
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( 28 ) 



In this equation, is the Neumann factor and F is the gamma function. Another 
difference between the FMM and the Anderson (1950) solution is that the number of 
modes associated with (p, denoted by m, is allowed to be greater than zero. 

The formulation incorporates a general spherical coordinate system. 
However, this coordinate system must be modified in order to accommodate different and 
more complex scatterer shapes. Conformal mapping is used to transform the Helmholtz 
equation from one set of coordinates to another in order to accommodate a variety of 
geometries. This modification is outlined below. 

b. Coordinate System 

The FMM uses an orthogonal coordinate system that can be created for a 
three-dimensional body of revolution from a two dimensional conformal mapping. The 
two dimensional approach was used by DiPema and Stanton (1994). The orthogonal 
coordinate system starts from the typical spherical coordinate system shown in Figure 2. 

The azimuthal angular coordinate is cp in this case, which ranges from 0 
to 2n, and is measured from the positive x-axis in the xy-plane. The polar angular 
coordinate is z?, which ranges from 0 to n, is measured from the positive z-axis. 


X 



Figure 2. Geometry for an irregular, axisymmetric, finite-length body of 
revolution. In this case, the body is symmetric about the z-axis. This 
figure shows the geometry for both spherical coordinates and an 
orthogonal, conformally mapped coordinate system, (From Reeder and 
Stanton, 2004) 
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The radial coordinate, r, is a scalar measured from the origin, which ranges from 0 to 0 °. 
The body of revolution is formed by rotating a contour about the z-axis in the (p 
direction. 

The new, orthogonal coordinate system has azimuthal angular coordinate 
V, which corresponds to (p and ranges from 0 to In. The polar angular coordinate is 
denoted by w, which corresponds to r? and spans a range from 0 to n. In the new 
coordinate system, the radial coordinate that defined the scatterer surface in the old 
coordinate system is changed. It is now defined by the locus of all points for which the 
new radial coordinate is a constant, or where u = 0 (Reeder and Stanton, 2004). The 
original body of revolution must also be parameterized in the new coordinate system. 
Using functions f{u,w) and g{u,w) together with trigonometry, the dimensions of the 
scatterer in the x, y and z directions become: 

x(u,w, v) = f(u,w) cos(v) 

(29) y{u, w,v) = f{u, w) sin(v) (Reeder and Stanton, 2004). 

z{u,w,v) = g{u,w) 

Further, the radial position vector becomes: 

(30) r{u,w,v) = x{u,w,v)i +y{u,w,v) j + z{u,w,v)k (Reeder and Stanton, 2004), 

where i , j and k are unit vectors in the respective x, y and z directions. The scatterer 
dimensions from Equation (29) can be substituted into this radial position vector to yield: 

(31) r = f(u,w)cos(v)i+f(u,w)sm(v)j + g(u,w)k (Reeder and Stanton, 2004). 

The orthogonal coordinate system is important because it aids in the calculation of 
normal particle velocities on the boundary, which are necessary in satisfying the 
boundary conditions (Reeder and Stanton, 2004). Orthogonality requires that r„ • r^, = 0, 

U, • T, = 0, and r ,• r, , = 0 . These conditions allow a conformal transformation, where a 
shape that was represented in the (x,y,z) coordinate system can be depicted as a shape in 
the {u,w,v) coordinate system. 
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c. Conformal Mapping and Parameterization 

In order to use the orthogonal eoordinate system, the {x,y,z) representation 
of the seatterer’s surfaee must be transformed with a mapping funetion into the new 
geometry in the {u,w,v) eoordinate system. Conformally mapping the old geometry into 
the new eoordinate system guarantees that the new geometry will be mutually orthogonal 
(Reeder and Stanton, 2004). Reeder and Stanton (2004) wrote “the eonformal mapping 
generates a new set of angular funetions whieh fit the seatterer surfaee more naturally; 
that is, points along the surfaee that ehange rapidly in (x,y,z) are plotted at a higher spatial 
rate yet are equally spaeed in (m,w,v).” 

The mapping funetion provides the means by whieh the seatterer shape is 
mapped to the new orthogonal, axisymmetrie eoordinate system. A three-dimensional 
mapping seheme does not yet exist within the field of mathematies, so the mapping must 
also be for an axisymmetrie body of revolution (Reeder and Stanton, 2004). To this 
point, the FMM has been developed with the two-dimensional mapping approaeh that 
was introdueed by DiPema and Stanton (1994). See Reeder and Stanton (2004) for 
further information and development of the DiPerna and Stanton (1994) mapping 
funetion. 

The mapping algorithm was not utilized in this researeh. For a smooth 
sphere or spheroid, the eonformal mapping reduees to: 

f = a*sin(w) 

(32) ^ 

g = b* oos(w) 

where a is the half-width of the body along the semi-minor axis, and b is the half-length 
of the body along the semi-major axis of the spheroid. Aspeet ratio {AR) is defined as: 

(33) AR = -. 

a 

The shapes that were used in eonjunetion with this researeh are shown in Figure 3. 
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Figures. A sphere with aspect ratio of 1:1 and a prolate spheroid with aspect 
ratio 5:1. These are the objects that are used in FMM simulations for 
this research. The objects are three-dimensional and are axisymmetric 
about the z-axis. 


With the change to the new coordinate system, the following relationships 
are established: 


(34) 


(p = v, 

r{u,w) = ^jf^{u,w) + g^{u,w), 


cos{'iKu,w)) = ^Y^. 

r{u,w) 


In addition, the Helmholtz equation becomes: 

(35) V^p{u,w,v) + k^F{M,w)p{M,w,v) = 0, 

where (u,v,w) make up the new coordinate system, and F{u,w) is a special function that 
depends upon the type of transformation (Reeder and Stanton, 2004). With this change, 
wave number becomes a function of position, but the Helmholtz equation itself is 
otherwise very similar to its counterpart in Cartesian coordinates (Reeder and Stanton, 
2004). 
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The solution to the new Helmholtz equation is: 


p“(.u,w,v)= I; X 

n=-o.m=-oo r{u,W) 

I 

n=-oom=-o<, r{u,W) 


(Reeder and Stanton, 2004). 


This equation yields the total far-field aeoustie pressure as a sum of the ineident and 
seattered aeoustie pressure in the new orthogonal coordinate system. This solution is 
valid for axisymmetric bodies for all frequencies and for all angles. It also applies to 
both impenetrable (soft and rigid) and penetrable (fluid) boundary conditions. 


d. Target Strength (TS) and Reduced Target Strength (RTS) 

At a great distance from the scatterer, the scattered acoustic pressure is 
given by the limit: 


(37) 



(Reeder and Stanton, 2004), 


where is the scattered pressure given by the second term of Equation (36), and p'"'" is 
the incident pressure given by the first term. The scattering amplitude, fs, is an 
expression of scattering efficiency and is a function of the scatterer’s size, shape, 
orientation, physical composition, in addition to the wavelength of the incident wave 
(Reeder and Stanton, 2004). The scattering amplitude is given by: 

(38) /. = Z Z i e'"- (Reeder and Stanton, 2004). 

„=-oom=-^ \r{u,w)) 


Acoustic energy that is calculated in the far-field region in the backscatter 
direction is usually conveyed as target strength {TS) in units of decibels (dB) relative to 
1 m (Urick, 1983). The expression for target strength is: 

(39) ra = 101oga-,,. 


where cr^HS the differential backscattering cross section. This is similar to the more 
common backscattering cross section, <7, but is multiplied by a factor of 4;rso that 
<7 = 4;r(7^^ (Reeder and Stanton, 2004). The differential backscattering cross section can 
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also be represented in the form cr^^=|/fc| , where is the scattering amplitude 
calculated in the backscattering direction (Reeder and Stanton, 2004). 

When comparing scattering events from similar objects of different sizes, 
it is often useful to normalize target strength by the square of an associated dimension. 
The resulting normalized target strength is then called “reduced” target strength (RTS). 
Using the length (L) of an elongated scatterer as the normalization constant, RTS can be 
expressed as: 


(40) 


RTS = \0{og 


f ry \ 

^bs 


\ L^ J 
1 


= 101og|/,,f-101og(Z^) 


= 101og 




bs 


L 


(Reeder and Stanton, 2004). 


For a sphere, the target strength should be normalized by instead of ll (Reeder and 
Stanton, 2004). Plots of RTS vs. dimensionless frequency {ka) are useful in resolving the 
scattering behavior of a scatterer, in addition to determining model characteristics and 
performance. RTS is expressed in units of decibels (dB). 

e. Benefits and Limitations of the Method 

Reeder and Stanton (2004) showed that the FMM is accurate over a wide 
range of scatterer cross sections, signal frequencies, and boundary conditions. Under 
many circumstances, the FMM is superior to other acoustic models. Figure 4 shows that 
the FMM matches the exact solutions for soft, rigid, and fluid spheres formulated by 
Anderson (1950). 

While the FMM has several advantages like enabling simulations of 
scattering from complex shapes, it also has several limitations. In order for the FMM to 
be applied to a three-dimensional scatterer, the scatterer itself must be axisymmetric, or 
symmetric about one of its axes. In other words, the outer boundary of the scatterer must 
be depicted by a function rotated around the length-wise axis. This is because the FMM 
was initially formulated by DiPerna and Stanton (1994) for the case of an infinitely-long 
cylinder, which is a two-dimensional scattering solution. Real scattering events typically 
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occur with three-dimensional seatterers. The FMM was developed from its original form 
by Reeder and Stanton (2004) to apply to three-dimensional bodies, with the limitation 
that the body must be deseribed by a function rotated about the length-wise z-axis. 





ka 


Figure 4, FMM comparisons to exact solutions. This shows reduced target 
strength in dB plotted as a function of ka for the cases of soft, rigid and 
fluid spheres. The exact solutions shown here are from Anderson 
(1950). (From Reeder and Stanton, 2004) 
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In addition, it is single-valued in r —there ean be only one value of r for eaeh value w 
(Reeder and Stanton, 2004). The FMM is also eomputationally expensive, in that it takes 
a long time to run on eurrent day eomputer systems. 

In addition to the above listed limitations, the FMM (or any seattering 
model, for that matter) is only as good as the mathematies that it ineorporates. The FMM 
uses spherieal wave funetions; however, not all of the seatterers to whieh this model is 
applied are actually spherical. For example, most fishes’ swim bladders are elongated, 
irregularly-shaped, prolate spheroids. Although the conformal mapping that is 
incorporated into the Reeder and Stanton (2004) FMM attempts to provide the details of a 
scatterer’s boundary, the basis functions within the formulation must still be somewhat 
similar to the object that they are being used to describe. Figure 5 shows the anatomy of 
an alewife (Alosa pseudoharengus) with the outline of the fish’s swim bladder from an x- 
ray image (Reeder, et al., 2004). This fish, and specifically its swim bladder and head, is 
an example of an individual scatterer to which the FMM would directly apply. 



Figure 5. X-ray image of an alewife {Alosa pseudoharengus) with swim bladder 
encircled in white. (After Reeder, et al., 2004) 


Notice that the swim bladder is not a sphere or a spheroid but a complex shape with an 
irregular boundary. Although the swim bladder is not spherical or spheroidal, the shape 
can be better represented by a prolate spheroid than with a sphere. The FMM would be 
more applicable to scattering from fish and submarines if it was formulated with 
spheroidal wave functions rather than spherical wave functions. 
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3. Boundary Conditions 

Inhomogeneities in the oeean may be characterized as hard-, soft-, or fluid-type 
scatterers, depending upon the material properties of which they are composed. A 
scatterer’s characterization depends upon the boundary between it and the ocean, as well 
as its interior structure. Sound interacts differently with each type of boundary and 
scattering body. Thus, there is a need for boundary conditions within a mathematical 
scattering solution, and boundary conditions often do not accurately depict the boundary 
of a scatterer. Different behavior is forced by each boundary condition, as shown by 
Figure 4 in the last section. It is important to use the correct boundary conditions, or the 
entire scattering formulation could yield incorrect results. 

Boundary conditions are used to produce the particular solution for a given set of 
initial conditions, such as scatterer properties. While Reeder and Stanton (2004) 
evaluated the model with soft, rigid, and fluid boundary conditions for their work, soft 
boundary conditions were used in this research for purposes of consistency in evaluating 
the effects of precision and numerical techniques. 

The soft boundary conditions are otherwise known as pressure release or 
Derichlet boundary conditions, where the total pressure vanishes along the surface of the 
scattering body (Reeder and Stanton, 2004). In other words: 

(41) p^\u„w,v) = Q 

from Equation (36). Thus, the series solution in the EMM general solution above is set 
equal to zero at the boundary. Reeder and Stanton (2004) outline the development of the 
boundary conditions for the specific case of a soft boundary, which is used here. The 
series coefficients, B , for the scattered field are: 

( 42 ) 

where the inverse notation indicates a matrix inversion and R and Q are integral 
expressions for summation over w (Reeder and Stanton, 2004). This boundary condition 
represents the physical case of a bubble or fish’s swim bladder, in which the difference in 
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fluid properties between the exterior and interior of the scatterer is signifieant and the 
surface of the scatterer is free to deform when insonified by an acoustic pressure wave. 

4, Series Solutions and Convergence 

For the FMM, the m-modes are the series components that are associated with the 
azimuthal angular coordinate (p. The n-modes are the series components associated with 
the polar angular coordinate z?. In a series solution, the modes are summed until a 
converged answer is achieved. The individual modes can be analyzed for their separate 
contributions to the total answer. When modes are added together, the solution should 
become increasingly converged with the addition of each mode. 

For the general formulation and solution to the spherical geometry in the 
Anderson (1950) paper, the series solution reaches convergence. The sum of the series 
components is an exact mathematical solution. For simple geometries, the series solution 
requires only a few modes in the summation to yield a converged solution. However, 
when irregular scatterer shapes are introduced into the problem, additional modes are 
required to be calculated in order to approach a converged solution. With more complex 
shapes, more modes are usually required to represent the scattering phenomena. 
Scattering from objects with higher aspect ratios or rough outer surfaces is extremely 
difficult to predict because of the higher number of modal combinations necessary in the 
computation. Based upon these circumstances, there is an inverse relationship between 
computational accuracy and computational efficiency. 

Convergence is achieved only to a certain point in the solution, however. When 
other factors such as error are introduced into the mathematical formulation by the higher 
modal combinations, model calculations become much more cumbersome. At higher 
modal combinations, error becomes a significant part of the solution. Error starts to 
dominate some modes for which the contribution to the total solution is only an 
extremely small number. In this circumstance, the error may outweigh the actual modal 
contribution, and the model will be inaccurate. At this point, the solution will incorporate 
numerous erroneous values that are not actually part of the intended solution, and the 
model can no longer be used for a useful scattering simulation. 
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5. Truncation 

Realistically speaking, a computed series can never give an exact solution because 
there are an infinite and unattainable number of terms in the solution. At a certain point 
within an infinite series, the series term being computed requires more decimal places 
than the computer can actually handle due to precision limits. The smaller terms at 
higher modes cannot be resolved by the computer. 

One method to deal with this problem is to truncate the solution and include only 
the modes that are known to be correct in the calculated solution. There is not a common 
stopping point in the series solution for all scatterer shapes. The point at which the 
solution is stopped, or truncated, is different for various initial conditions and scattering 
surfaces. 

Truncation may also be required when the time needed to calculate the next mode 
would otherwise extend model runtime beyond the useful limit. In many situations, 
truncation at a certain point in the series solution could be necessary, in order to avoid 
computer delays for additional calculations that yield little or no benefit. There is little 
value added in calculating a mode that changes the solution only by a fraction of a 
percentage point. 

Reeder and Stanton (2004) conducted a study to discern the operational range of 
the FMM. They tallied the results in the performance envelope plot shown in Figure 6. 
They defined a converged solution as “one in which the computation of additional modes 
does not significantly change the result for a given value of ka" (Reeder and Stanton, 
2004). The truncated approximations are those for which the FMM yields results that 
leave some question of reliability and accuracy. The numerically stable approximations 
are those that reach beyond the truncated approximations and that still yield results that 
are not rendered completely useless by computational errors. 
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Soft Spheroid Performance Envelope-Broadside Backscatter 



Aspect Ratio 

Figure 6. Performance envelope for a broadside incident plane wave for a smooth 
prolate spheroid under soft boundary conditions. The performance 
level is indicated by the shading and is plotted as a function of ka and 
aspect ratio. (From Reeder and Stanton, 2004) 

Numerical techniques are implemented in this research with two principle goals: 
(1) extending the aforementioned performance envelope, and (2) making the lower modal 
combinations more accurate and, therefore, more useful. If either of these goals is 
achieved, then model performance is improved. The intention is to investigate changes in 
the performance envelope based on targeted changes in the formulation or calculation of 
results. 

6. Accuracy versus Precision when Applied to Computers 

In the course of this research, only a sphere and a prolate spheroid were 
considered as scattering shapes. For these shapes, the accuracy of the FMM was 
previously demonstrated by Reeder and Stanton (2004) in comparison to the Anderson 
(1950) exact solution, as shown in Figure 3, and other models for a prolate spheroid. 
While Reeder and Stanton (2004) addressed the problem of accuracy, they did not 
entirely address the issue of precision. Precision is inherent in any formulated computer 
calculation. Answers do not change from one model run to another, unless the inputs to 
that model are altered. However, another form of precision is considered here, which is 
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known as computer precision. With different values of precision, changes in the numbers 
of binary digits and decimal digits utilized in calculations will craft different answers. 

The precision of a particular computer is only one facet of the computing process. 
The following sections reveal the complexity of computer calculations, while introducing 
the key computer fundamentals behind this research. 

D. APPLIED COMPUTER TERMS AND PROCESSES 

Nearly all numerical computing utilizes floating-point arithmetic. Almost all 
computers use the Institute for Electrical and Electronics Engineers (IEEE) binary 
floating-point standard to represent numbers for computer storage. In the late 1940s and 
early 1950s, many people were afraid that the intrinsic rounding errors of floating-point 
computing would make the results of complex calculations too inaccurate to be useful for 
science (Overton, 2001). That is still the fear for the EMM at higher frequencies and 
modes for complex shapes. In order to further understand this problem and to elucidate 
the scope of this thesis, some introductory material in floating-point, computer precision 
and other computer terminology is included. 

1. Introduction to Floating-Point Arithmetic 

A “bit” of computer storage space is a single binary digit. A “byte” is a group of 
8 bits. A “word” is 4 consecutive bytes of computer storage space, or 32 bits, while a 
double-word is 8 consecutive bytes, or 64 bits. All real numbers have binary 
representations. 

Eloating-point representation is based primarily upon scientific notation. In 
floating-point form, the number S is called the significand while the number E is called 
the exponent. Eor example, the decimal number 0.00000625 can be expressed in 
scientific notation as 6.25 x 10'^ in base 10 notation, where 6.25 is the significand and 10' 
^ is the exponent. This is called floating-point because the decimal point in the decimal 
number “floats” to the position immediately after the first nonzero digit in the decimal 
expansion of the number (Overton, 2001). However, the computer standard is base 2, so 
the above expression would be written in the form x = ±5* x 2^, where 1 < iS < 2, and the 
binary point floats to enable representation in this form. This representation is called the 
normalized form. In order to store a floating-point number, the computer word is divided 
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into three values, representing the sign of the number (i.e., 0 for a positive number and 1 
for a negative number), the signifieand S, and the exponent E. For example, a 32-bit 
word could use 1 bit for the sign, 23 bits for the signifieand, and 8 bits for the exponent 
of a floating-point number (Overton, 2001). The binary expansion of the signifieand 
would then be: 

(43) S = {b,.bp,b,..b,,) 

where b^ is always 1. 

It is important to note that the bit on the left of the binary point (hg) will always 

have the value 1 in binary representation for a floating-point number. In order to avoid 
using valuable storage space, this value is not stored and it is called the hidden bit 
(Overton, 2001). This is an important gain in precision computing, because it allows the 
storage of another binary digit. 

While all real numbers have a binary representation, not all of them can be stored 
in a prescribed floating-point form. For a given amount of storage space, not all numbers 
can be represented within that numerical capacity. For the 32-bit word mentioned above, 
the exponent must fit in the range -128 < £' < 127 because all 8 bits of storage space for 
the exponent would be filled with the binary representation of E, which is 01111111 for 
E=\21 (Overton, 2001). A number larger than 127 would require more than 32 bits of 
storage space, which is not allowed by 32-bit computer architecture with only 8 bits 
prescribed for the exponent. Many numbers must be rounded before they can fit into 
floating-point form because the binary expansion must contain an exponential power of 2 
within a prescribed range. This rounding process introduces error, which will be 
discussed later. Even the number 1/10 does not have a finite binary representation and 
will introduce error into a computing process (Overton, 2001). 

2. Precision and Machine Epsilon 

Within a given bit string, numbers can be prescribed for the signifieand and 
exponent. The number of bits stipulated for the signifieand (plus one for the hidden bit) 
is called the precision of the floating point-number. So, for the system described above 
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where S=23, the preeision is 24. For any floating-point number x with precision p, 

(44) x = ±(l.hA...V,V,),x2". 

This means that the smallest floating-point number greater than 1 is then: 

(45) (1.00...01)2 =l + 2-<^-'). 

The difference between this number the number 1 is called machine epsilon, which is 
usually denoted as eps, or g . This can be expressed as: 

(46) G= (0.00... 01)2 =2-(^-'^ 

Machine eps is the smallest number that can be discerned by a computer, which is also 
the smallest difference between numbers that can exist on any given machine. With each 
different value of machine precision, eps changes. So, by changing machine precision, 
which changes the number of binary digits in the storage space of S, the machine eps is 
also specified. However, machine precision cannot be changed on a given computer 
unless different hardware is installed to facilitate changes in numerical storage 
capabilities. 

Higher values of precision allow a larger range of numbers to be stored within a 
computer. The largest number that can be stored by a computer is increased, while the 
smallest number that the computer can store is decreased. Higher precision therefore 
leads to increased technical capabilities for a computer system. 

Precision represents a challenge for certain calculations. Most computations 
require only a basic level of precision; however, certain computations require a large 
degree of precision in order to yield results that are accurate enough to be useful. It is 
believed that the FMM is currently limited by machine precision for higher modes and 
for shapes of high complexity and eccentricity. According to Reeder and Stanton (2004), 

The FMM generates a transition matrix, much like the T-matrix model 
(Waterman, 1968), that relates the incident field coefficients to the 
scattered field coefficients. For a spherical scatterer, the transition matrix 
is diagonal and each nonzero term on the main diagonal is an eigenvalue 
for each mode computed. If the scatterer shape deviates from spherical, 
the matrix contains off-diagonal terms. The additional higher modal terms 
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required to represent the scattering become extremely small, sometimes 
falling below the value that can be accurately represented numerically, 
resulting in a singular matrix in which the true values of its elements are 
below the precision of the machine. (Reeder and Stanton, 2004) 

The resurgent problem is that the actual elements of the ill-conditioned transition matrix 
in the FMM require more binary digits and smaller machine epsilon than 32-bit precision, 
and even 64-bit precision, will allow. 

Before proceeding, it is important to note that there are two forms of computer 
precision—machine precision and software precision. Machine precision is a set value 
on any given computer, based on how that machine was built, and the value of precision 
that is most commonly referenced in literature. Machine precision governs how fast a 
computer will perform a calculation, based on the number of bits used for storage. 
However, certain software programs introduce another facet to the concept of precision. 
Tools have been developed to allow software precision changes on the fly that will allow 
certain computations to be carried out at lower or higher precision than the machine 
architecture would otherwise allow. Software precision is conceptually identical to 
machine precision but can be changed with clever computer coding. MATLAB is an 
example of a software program that can run calculations at a higher precision than that of 
the machine on which it runs (i.e., MATLAB can run calculations at double precision, 
even if it is running on a computer that was built with single precision). For example, the 
version of MATLAB employed in this research uses two 32-bit words to store a single 
floating-point number, effectively gaining 64-bit precision via software syntax that 
instructs the computer to do so. 

3. Single, Double, Extended, and Variable Precision 
The IEEE standard has two fundamental formats, which they call single and 
double. The 32-bit word is the typical mode of storage for computers. Single format 
numbers use a 32-bit word as their storage mechanism. However, single format does not 
offer enough binary digits for applications where higher precision is needed or where a 
greater exponent is required. Double format uses a 64-bit word, which offers a larger 
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number of binary digits for the significand and exponent of the normalized 
representation. For double format, the exponent may take values within the range: 

(47) -1022 <.£<1023. 

Table 1 displays a summary of the ranges of E and the maximum (A^max) and minimum 
(A^min) numbers that single and double formats ean handle. Table 2 shows how double 
precision numbers are stored within the exponent and significand. 


Format 

E^min 

Emsa 

Amin 

Amax 

Single 

-126 

111 

2-'^^ = 1.2 X 10'^^ 

= 2'^^ = 3.4 X 10^^ 

Double 

-1022 

1023 

2'^°^^ = 2.2 X 10'^°^ 

= 2^°^^ = 1.8 X 10^°^ 


Table 1. Range of Numerical Values for IEEE Floating Point Single and Double 
Precision Formats. (From Overton, 2001) 

The IEEE standard also has an extended format, which holds 15 bits for the 
exponent, 63 bits for the significand, one digit for sign, and one for the leading (hidden) 
bit that is not hidden in this format (Overton, 2001). Machines that run the extended 
format often run more slowly than their single and double precision counterparts, because 
every numerical value is stored with extended precision. Table 3 shows the precision and 
eps values for the various IEEE standard floating-point formats. 

Machine precision cannot be changed without changing the actual hardware 
storage devices that are built into a computer. In order to circumvent the current limits 
on machine precision, symbolic representations of variables and variable-precision 
arithmetic (VPA) can be used to increase accuracy and precision for certain calculations. 
Symbolic mathematics that enable variable precision and rational representations of 
numbers are discussed in more detail in Chapter III. 
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± 


C11C12C13...C11I 


b]b 2 b 3 ...b 52 


If exponent bitstring is ai...aii 

Then numerical value represented is 

( 00000000000 ) 2 = (0)io 

( 00000000001 ) 2 = (l)io 
( 00000000010 ) 2 = (2)io 

(00000000011)2= (3)io 

i 

(01111111111)2=(1023)io 

± (0.hih2h3...h52)2 X 2'^°^^ 

± (I.hih2h3...h52)2 X 2'^°^^ 

± (I.hih2h3...h52)2 X 2'^°^^ 

± (I.hih2h3...h52)2 X 2'^°^“ 

i 

± (I.hih2h3...h52)2 X 2° 

(10000000000)2= (1024)io 

i 

(11111111100)2= (2044)io 

± (I.hih2h3...h52)2 X 2 ' 

i 

± (1 .bibib^.. .bsi)! ^ 

(11111111101)2= (2045)io 

± (1 .bib2b3. ..^ 52)2 X 2 

(11111111110)2= (2046)io 

± (I.hih2h3...h52)2 X 2 '“^ 

(11111111111)2=(2047)io 

± 00 if hi = ... = ^52 = 0, NaN otherwise 


Table 2. IEEE Double Format with E represented by and S represented 

by bib 2 b 3 ...b 52 . The left-hand column shows binary and decimal 
representation of the exponent strings. The right-hand column shows the 
normalized representation of the numerical value. (From Overton, 2001) 


Format 

Precision 

Machine Epsilon 

Single 

p=24 

G =2'^^ = 1.2 X lO'^ 

Double 

cn 

II 

G = 2'^^ ~ 2.2 X IQ-'® 

Extended 

II 

G =2'“ ~ 1.1 X lO-'" 


Table 3. Precision of IEEE Floating-Point Representations. The precision, p, is 
the number of binary digits in the significand of the normalized form of 
the stored value. Machine epsilon is the smallest distinguishable number 
associated with the given format. Although the precision seems high in 
binary form, the decimal equivalent of these numbers of binary digits of 
precision is much lower. For example, double precision floating-point 
numbers in binary representation with 53 bits of precision must be 
rounded to their decimal equivalents with only 15 or 16 decimal digits. 
Thus increasing the number of decimal digits used in computation can 
make a much greater impact than increasing the number of binary digits 
by the same value. (From Overton, 2001) 
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4. Computational Expense and Realistic Limits 

The term “expense” refers to the amount of eomputational time and computer 
storage that is required for computer operations. Computational storage increases with 
higher precision, because more bits are required for storage of larger bit strings. 
Computational time requirements also increase as precision increases, since more I’s and 
O’s need to be placed into the bit strings with each operation. At a given value of 
precision, computational time can only be lessened with the use of a faster computer 
processor or a more efficient mathematical formulation. This means that the potential use 
of higher precision is reliant upon the amount of time available for computer operations. 
Higher precision is more expensive to the end user that is cognizant of limited resources. 
E. COMPUTATIONAL ERROR 

1. Introduction to Computer Error 

As mentioned above, rounding takes place when a number cannot be stored 
within a bit string of prescribed length. Numbers are rounded when converted from 
decimal to binary for computation, and they are again rounded when they are converted 
from the binary storage form to decimal screen output. In addition, when numbers are 
added, subtracted, multiplied or divided, the result might be a number that cannot be 
represented in floating-point without being rounded. This is called integer overflow, 
which the computer corrects by rounding the number before it can be stored. This 
compounds the error that may have started with numbers that were not floating-point 
numbers to begin with and required rounding. 

In a simple addition such as: 

1/2 + 1/3, 

there are really three roundoff errors associated with the output. Note that 1/2 is actually 
a power of 2 and does not require rounding prior to the binary storage of the number. 
The first roundoff error is encountered with the division of 1 by 3 when the computer has 
to round in order to store the result of the division within the prescribed number of bits in 
the significand of the associated IEEE format. Another error is introduced with the 
addition of 1/2 to the result of the aforementioned division, because this result must also 
be rounded and stored. The final error occurs when the binary result is converted to 
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decimal format for the printed output. When IEEE double floating-point arithmetic is 
used to calculate this result, 64 bits are used for the stored result, but only 53 are used in 
the binary expansion of the significand. This usually gives about 16 decimal digits, when 
the conversion is made. However, only 15 digits print to screen in this specific case, with 
the result showing: 

0.83333333333333 (The Mathworks, Inc., 2002). 

This result is not the exact answer to the problem above because it incorporates the three 
errors that were introduced during its creation. 

2. Compounding Error - The Importance of Precision 
As already mentioned, a calculation performed below a required level of precision 
will incorporate roundoff error into the answer. This occurs because the real values have 
more decimal or binary digits than the computer can store within its floating-point 
representation. The mathematical solution becomes limited by machine precision, so that 
the computer solution does not accurately represent the actual intended computation. If 
this rounded answer is then used in further calculation, it can propagate and multiply the 
roundoff error by performing various mathematical operations on an already erroneous 
number and heavily impact the accuracy of the final solution. This happens frequently in 
models such as the EMM, in which the smallest of numbers is important to the 
mathematical computation of an accurate answer. 

Small changes in initial conditions within a model should produce small changes 
within the results. In this situation, the model is said to be stable. However, if a small 
change in initial conditions produces a drastically different answer, then the model is said 
to be unstable (Mathews and Eink, 1999). The errors discussed above may introduce 
instabilities throughout the course of model calculations. 

F. NUMERICAL TECHNIQUES 

1. Operator Modifications and Reducing Roundoff Error 
Addition, subtraction, multiplication and division all inherently incorporate 
different amounts of roundoff error. With each calculation, there is an associated error 
either for storage or output to the computer screen. An efficient mathematical 
formulation reduces the number of operations that must be performed. When fewer 
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mathematical operations are performed, fewer incidences of rounding occur which 
consequently minimizes roundoff error. 

2. Singular Value Decomposition 

Singular value decomposition (SVD) is a mathematical tool that assists in 
handling series of equations or matrices that are singular or nearly singular (Press et al, 
1992). According to Press et al, “Any M x N matrix A whose number of rows M is 
greater than or equal to its number of column s N, can be written as the product of an M x 
N column-orthogonal matrix U, an A x A diagonal matrix W with positive or zero 
elements (the singular values), and the transpose of an Ax A orthogonal matrix V.” This 
technique is extremely useful in the case of an ill-conditioned matrix that includes 
elements below machine precision. The SVD technique sets elements in a matrix below 
a certain threshold to zero. 

Essential to the EMM are several values within the transition matrices that are 
extremely small. These values can make important contributions to the solution if they 
are real parts of the solution, but they are often introduced into the problem by roundoff 
error. In the case of the Reeder and Stanton (2004) model, SVD removes all elements 
below a manually set threshold that are zero or near zero within the calculated matrices. 
“Singular values whose ratio to the largest singular value is less than A times the machine 
precision are set to zero” (Press et al., 1992). Problems exist when this threshold is set 
too high or too low. If it is set too high, then some values will be removed from the 
matrices that are actually valuable parts of the scattering solution. If the threshold is set 
too low, many values that are supposed to be zero but are non-zero because of roundoff 
errors may be left in the solution as erroneous contributions to the scattering solution. 
The threshold for SVD should be set at-or-below the value of eps associated with the 
increased precision. This will prevent the gains of increased precision from being 
removed by the SVD commands, which conduct valuable quality control for model 
calculations. 

Singular value decomposition can be a valuable and powerful tool when used 
properly. Its intention is to remove “erroneous subspaces” in order to produce a more 
stable result (Reeder and Stanton, 2004). If used incorrectly, SVD can eliminate matrix 
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values that are positive contributions to the solution. Incorrect implementation could 
result in reduction of amplitude and damage to the composition of the solution. 

G. EXPECTED RESULTS 

The planned model improvements should reduce or eliminate model instability 
that occurs in many of the FMM’s higher modal combinations due to compounding 
roundoff error. The roundoff errors are inherent in nearly all floating-point numbers in 
the model calculations, and the improvements made in this research are aimed at 
decreasing the roundoff error to a minimal amount. The improved model should have 
less instability in the higher modal combinations, in the higher frequency ranges, and 
when applied to shapes of greater eccentricity and complexity. 

Increased accuracy is another anticipated attribute of the improved FMM. If the 
increased precision can add extra decimal places or binary digits to the values used in 
calculations, then the expectation is to see greater accuracy in the results of the improved 
model in comparison to the results of the FMM used in Reeder and Stanton (2004). 

Implemented techniques will also aim to improve the performance envelope of 
the FMM. If the model is more accurate, fewer modes will be required to arrive at a 
converged solution. This may reduce the overall number of modes that are needed to 
gain useful results. 
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III. IMPLEMENTATION OF THEORETICAL IMPROVEMENTS 


The following sections describe the background work that was required to 
implement various numerical techniques, including information on hardware, software 
and coding routines that were used to gain numerical results. This chapter covers the 
specific tools within MATLAB and the numerical techniques that were used to effect 
changes in mathematical calculations in the FMM, in addition to the methods for 
evaluating the research results. 

A. TERMINOLOGY AND RESEARCH SETUP 

In the following discussion of implementation and results, “original” model refers 
to the model developed by Reeder and Stanton (2004). Additionally, the “improved” 
model refers to a specific version of the FMM that incorporates the individual changes 
that are explicated in the following sections. In most cases, the improved model is 
compared to the original model through direct comparison of output plots, in order to 
ascertain the added value of the executed technique. The exact formulations to which 
Reeder and Stanton compared their results, such as the Anderson (1950) solution, are not 
used in this research, because they too incorporate roundoff errors when calculated at 
double precision. 

B. HARDWARE AND SOFTWARE SPECIFICATIONS 

1. Hardware 

The machines used for the bulk of this research are Dell Precision 670 n-series 
machines with Xeon processors. These machines are built to use single precision with 32 
bits of storage space for each floating point number. Although these are not the fastest or 
most robust machines on the market today, they do represent the average technical 
computing machine that would be used in a non-research, or operational, environment. 
This provides an element of realism to the research, since the results can be directly 
applied to the operational environment that uses similar machines. Most machines that 
are currently being built incorporate 64-bit machine precision, so the research platforms 
are already surpassed by the storage capabilities of the newest systems. 
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2. Software 

All computer code generated for this research was created or modified in 
MATLAB Version 7.0.1.24704 (R14) Service Pack 1, which is produced by The 
Mathworks, Inc. The 3.1.1 (R14) Service Pack 1 version of the Symbolic Math Toolbox 
was used for all work associated with symbolic numbers and variable-precision 
arithmetic. 

Although the machines described above run at single precision, MATLAB stores 
variables and performs calculations in double precision, doubling the actual machine 
precision through the software syntax. This was verified using the eps command within 
MATLAB, which displayed eps values that corresponded to double precision on a 
machine that was known to run at single machine precision. 

The Dell machines used in this experiment ran Linux Kernel 2.4.21-32.0.l.ELsmp 
#1 SMP May 17 17:52:23 EDT 2005 1686 for all applied research. The Einux operating 
system is very supportive of technical computing and proved to be a reliable research 
tool. 

C. INCREASING PRECISION THROUGH SYMBOLIC MATHEMATICS 

The Symbolic Math Toolbox within MATLAB encompasses over 100 functions 
that allow exploitation of the Maple kernel through the use of specific syntax in the 
MATLAB language (The Mathworks, Inc., 2002). This enables symbolic calculations to 
be performed within the floating-point MATLAB arena. This is an extension of the 
capabilities of MATLAB that enhances the standard computing and plotting tools. 

The use of the functions within the Symbolic Math Toolbox creates a variable of a 
new class, called symbolic object, or sym for short, which is not a floating-point number. 
It is a symbolic representation similar to that created when doing calculations with pencil 
and paper. Internally, the computer stores the sym as a character string. Syms can be 
used to represent symbolic variables, expressions, and matrices. To obtain the numerical 
value of a sym, the double command must be used, which converts the sym into a 
floating-point number. It is important to note that this conversion requires the number to 
be rounded to the nearest floating-point value. MATLAB also performs calculations 
differently with values of class sym than it does with floating-point numbers of class 
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double. Operations performed with double precision will yield answers in double. 
Operations performed with syms will yield answers with rational form as syms. 
MATLAB actually searches for the lowest common denominator of the fractions 
represented by the sym numbers and performs calculations in rational arithmetic, with the 
help of Maple (The Mathworks, Inc., 2002). 

Three general representations of syms are possible within the Symbolic Math 
Toolbox utility, including numerical (regular floating-point), variable-precision 
arithmetic (user-defined), and symbolic representation (rational/exact). These 
representations are all different, but are all class sym variables. Floating-point, or 
numerical, operations are the least expensive in computational time and memory out of 
the three forms of sym variables, but the results are not exact due to associated roundoff 
errors. Rational operations use exact rational values without any error in their 
representation, given that the numbers required for calculation can be represented as 
ratios or integers. MATLAB stores sym numbers in rational form (i.e., 1/5), rather than 
using floating-point representation. However, symbolic numbers take the most computer 
time and memory to compute out of any of the available numerical forms in MATLAB 
(The Mathworks, Inc., 2002). Calculations performed with variable-precision arithmetic 
(VPA) are between the other two forms with regards to expense, but may afford 
somewhat of a balance of accuracy and computational time for some calculations. The 
following paragraphs describe the two specific symbolic forms that were used in this 
research. 

1. Variable-Precision Arithmetic 

One way in which machine precision can be extended is to use VPA. With the 
vpa command that is part of the Symbolic Math Toolbox in MATLAB, a user can set 
precision to a desired value for certain calculations. VPA uses the digits command to set 
the number of decimal digits of precision at which a calculation will be made with 
symbolic math and variables of class sym. Mathworks, Inc. says that the vpa command 
uses “variable precision floating-point arithmetic with D decimal digit accuracy, where D 
is the current setting of digits" (MATLAB Help Document for vpa.m). The resulting 
expression from vpa is a symbolic value or array that is similar to a character string. 
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With digits=\Q, the number 1 is represented by VPA as “1.000000000” throughout all 
subsequent ealeulations, with the value 1 eounted as one of the ten digits set by digits. 

The default value of digits, unless otherwise speeified, is 32 deeimal digits. In 
other words, the number “1” would be expressed as the number 1 followed by 31 zeroes 
after the deeimal point. This inereases the aeeuracy that ean be earned through 
ealeulations by avoiding roundoff error until the last decimal digit of a long string of 
decimal digits, whereas double precision introduces roundoff error to store as binary and 
then rounds again to convert to sixteen decimal digits of output. Thus, roundoff error is 
still an unavoidable part of VPA calculations, but it is effectively reduced for calculations 
in which fewer than D digits are required to accurately represent the solution. 

Even by using VPA with digits=\6, which is the same as the MATLAB default 
output in double precision, the outcome of certain calculations may be improved. This is 
because VPA values with increased numbers of decimal digits are stored and carried 
through calculations as character strings, rather than using binary storage and floating¬ 
point operations. However, VPA must be used in a way that exploits precision to gain 
accuracy in smaller decimal places. Otherwise, computational expense is added with no 
numerical benefit. In the EMM, higher modal combinations use this extended precision 
afforded by VPA. VPA is the primary tool used in this research. 

Implementation of VPA into the Reeder and Stanton (2004) EMM introduced a 
multitude of programming errors that had to be overcome. Commands such as max, find, 
and even sqrt can not handle the symbolic variables that are created with the vpa utility. 
Some of these difficulties were overcome via the creation of “dummy” variables within 
MATLAB. The dummy variables are mirror images of the original variables passed to 
the max, find, and sqrt commands. Before the original variables are passed to these 
commands in which errors occur, they are converted from symbolic form back to double 
precision floating-point format by using the double command. The variable can then be 
reset to its original precision after bypassing the command in which it would have 
otherwise produced an error. This is accomplished by setting the original variable name 
equal to the dummy variable. By using this technique, subsequent lines of code do not 

have to be modified to accommodate new variable names or classes. 
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The VPA techniques described above were the primary means of improving the 
FMM. The vpa command was used throughout the entire model, where possible, 
increasing the precision of the majority of calculations performed by the FMM. 

2. Symbolic Variables and Rational Arithmetic 

The sym command within MATLAB creates variables that can be of a few 
different forms: floating-point, rational, expressions in terms of eps, and decimal 
representations with the number of decimal digits set by digits. The rational numbers 
created by the sym command introduce an even higher level of accuracy into computer 
calculations at the expense of significantly greater computational time. 

In this research, the rational representation of sym, which is the default used by 
MATLAB if not otherwise specified, was used to generate exact rational numbers. In 
this way, calculations performed with sym are done with exact numbers and avoid 
roundoff errors inherent in floating-point and vpa calculations. If a number cannot be 
represented in rational form, then the result will be an irrational number in floating-point 
format with default precision. 

Symbolic math with rational numbers eliminates any and all roundoff errors 
associated with the symbolic variables, because the rational numbers do not require 
rounding. The sym command would produce exact results for the FMM, if all values 
could be carried through the model as rational numbers. The sym command was 
implemented into the FMM to evaluate another of the symbolic variable formats, but the 
computational expense that was inherent in the associated calculations was too great to 
compute even the first mode in the FMM solution. The model would not run to 
completion. Errors were encountered before the model could produce results even at 
mode m=\, n=l. These errors are described at the end of this chapter. The rational 
numbers produced by sym are known to be the most expensive of all variable forms 
within MATLAB, so the errors experienced here were not surprising. 

D. NUMERICAL TECHNIQUES 

1. Operator Modifications 

With each calculation that is conducted within MATLAB, there is inherent 
roundoff error. The simplest technique in reducing roundoff error is to reduce the 
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number of required calculations. In order to do this, all calculations that can be 
simplified are reduced to their most basic form. Reducing the total number of 
calculations, even if only by a few, reduces the number of rounding errors incorporated 
into and amplified within the solution. In certain calculations, the roundoff error could be 
adding a large enough amount of error that it could affect the solution quite significantly. 
The FMM code was carefully modified in this research to reduce the total number of 
calculations being performed by the model. 

2. Singular Value Decomposition 

Singular value decomposition was not treated in this research because the modal 
combinations where SVD made significant changes to the answer in the original FMM 
were outside the capable range of the improved FMM. The range of the improved FMM 
will be discussed more in Chapter IV. The SVD processes that were incorporated into 
the model by Reeder and Stanton (2004) remained the same for this research. 

E. EVALUATION TECHNIQUES 

1. Visual Inspection of Reduced Target Strength (RTS) Versus Non- 
dimensional Frequency (ka) 

Results from the original Reeder and Stanton (2004) model are plotted 
simultaneously with the results of the improved model. Graphical comparison of the 
solution achieved with the improved model to the results of the original model shows 
improvements in the solution. 

2. Measurements of Computational Time 

Computational time was recorded for each model run with the use of the tic and 
toe commands within MATLAB. The time was started at the beginning of each model 
run with tic, and the runtime was recorded at the conclusion of the calculation of each 
modal combination with toe. In this way, a cost-benefit analysis could be conducted on 
the results of the model runs, based upon the model’s runtime and performance gains. 

F. MATLAB/CODE LIMITATIONS 

One of the major limitations on the use of VP A was the symbolic form of the 
output. Calculations with symbolic numbers could be conducted to a certain extent 
before errors were encountered in MATLAB in association with the variables of class 
sym. The symbolic numbers had to be evaluated with the eval or double command in 
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order to convert them back to double precision so that they could be used in further 
calculations beyond the point of the errors. This is desirable because the higher precision 
values could be carried through most of the model calculations, and the more accurate 
results from those calculations could then be used in further computations in double 
precision. 

In addition to the max, find, and sqrt commands mentioned earlier, other 
operations and commands that posed challenges to the implementation of symbolic 
variables included indexing operations within matrices and saving variables in symbolic 
format. The squeeze command that is used within the integration process in the FMM 
would also not work with variables of class sym. This meant that the implementation of 
the VPA or sym variables could only be carried through the end of the calculations of the 
Bessel functions and associated Legendre functions. Once the matrices for these 
functions were populated by the model, the results had to be converted back into double 
precision for further computation with the squeeze command. The errors in this 
paragraph were encountered when sym variables were introduced into lines of computer 
code containing the commands listed above. This meant that the variables had to be 
changed to double precision to carry them through these lines of code. 

The most significant error was one that limited the number of modes that could be 
computed by the improved FMM once VPA was already implemented. An 
insurmountable error encountered within MATLAB at high modal combinations while 
running calculations with vpa variables was displayed as “Error, integer too large in 
context” in the MATLAB command window. This error corresponds to Mathworks 
Technical Solution number 1-1AG3M on the Mathworks.com support website. The 
solution reads: 

Our development staff has been notified and is currently looking into 
addressing this problem in a future release of the Symbolic Math Toolbox. 

This may be a problem with the way memory management is performed 
by Maple. At present, the only potential workaround is to wrap your call 
to a symbolic calculation that operates on numbers with large numbers of 
digits inside a TRY/CATCH block and clear the Maple function... (The 
Mathworks, Inc., 2005) 
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This statement at least shows promise for the next-generation Symbolic Math Toolbox. 
However, the “workaround” presented here does not solve the problem. It merely steps 
around the error message and proceeds without the numbers that generated the error in 
the first place. Technical improvements may soon be incorporated into MATLAB that 
would allow some of the roadblocks associated with the symbolic mathematics in this 
research to be avoided in future work with the FMM. 
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IV. NUMERICAL RESULTS 


A. GENERAL CHARACTERISTICS OF THE IMPROVED MODEL 

There are several speeifie changes and traits to note in the improved model, but 
the general changes are noted first for perspective. First, the improved FMM requires a 
much greater amount of computational time than the original FMM. This was expected 
because of the use of symbolic variables, in addition to the increased precision. Second, 
very few modes could be calculated with the improved model because of the 
computational errors inherent in MATLAB, which were discussed in Chapter III. The 
model encountered multiple errors, mostly due to several MATLAB functions’ 
incompatibility with the symbolic variables created by the vpa and sym commands. 
Third, the research exposed many findings during the implementation of VPA that have 
been left as items for follow-on research. The following paragraphs expound upon some 
of these general characteristics and offer insight into the model output that was obtained 
with the improved model. 

B. RESULTS OF EXTENDED PRECISION AND NUMERICAL 

TECHNIQUES 

A variety of techniques were implemented to improve the FMM. The applied 
operator modifications for addition, subtraction, multiplication and division operations 
did not produce any noticeable changes in the results of the original FMM. This is most 
likely because there were no apparent divisions by extremely small numbers or 
multiplications of extremely large numbers where the last few significant digits of the 
computed values make a large contribution to the end result. Although no improvements 
were observed in the results of the improved FMM based on the operator modifications 
alone, the code was modified for efficiency to ensure that only the smallest number of 
calculations required would be conducted by the improved FMM when precision was 
later increased. Increases in precision were the next step in improving the model. 

The results of the improved FMM are revealed in the following paragraphs and 
plots. Model updates used in producing these results include VPA and operator 
modifications. The effects of SVD are discussed as a point for further research. 
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Through several iterations of attempted VPA implementation, the TS vs. ka plots 
did not change. This information suggested three possible situations: (1) VPA was not 
being implemented correctly, (2) the improved results from the VPA implementation 
were being negated by roundoff or SVD, or (3) VPA did not make a difference in model 
performance. The first two of these situations were investigated in-full to guarantee the 
validity of the results. This led to several changes in the VPA implementation, which 
resulted in several discoveries in the way that vpa actually works with the computer- 
stored variables. The stored variables were checked after each line of code to guarantee 
that the specified level of precision was being carried through calculations. Once VPA 
was determined to have been implemented correctly, the model was run again, and it 
produced significant differences from the output of the original model. The VPA 
technique was the most effective of all of the executed numerical improvements. 

Figure 7 shows the first modal combination of the improved FMM with digits=\6, 
which produced the same numerical results as the original FMM. The fact that the 
improved model curve is exactly the same as the original model curve for the initial 
modal combination of m=\, n=\ gives confidence in the execution of VPA. 

The improved FMM yielded results for only the first two modal combinations 
(i.e., m=\, n=\ and m=\, n=2) with digits=\6 for a soft sphere. This gave an initial 
indication that the improved FMM would be extremely limited in the range of modes that 
it would compute for higher values of digits. The improved FMM was run with 
digits=24, digits=32, and digits=40, corresponding to 1 14, 2, and 2 14 times the precision 
of the original FMM, respectively. These runs yielded results for the first two modal 
combinations for all values of digits and results for the third mode (i.e., m=\, n=3) for 
digits=32 only. Because results for three modal combinations were attained for digits=32 
and because 32 digits corresponds to exactly twice the precision of the original FMM, the 
results for digits=32 were used for further analysis. 

The improved FMM results for digits=32 are shown in Figure 8. For m=\ and 
n=2, shown in Figure 8(b), the improved model deviates from the original. The shape of 
the RTS vs. ka plot for this modal combination suggests that the addition of one mode at 


50 



higher precision improves the result of the original model at the same modal combination 
by extending the converged, smooth portion of the curve to a higher value of frequency, 
or ka. The same holds true for mode m=\, n=3, which is displayed in Figure 8(c). The 
extension of the smooth portion of the curves in these plots indicates that the interference 
patterns in the representative mathematical functions are less destructive to the solution at 
lower values of ka. With less destructive interference at lower values of ka, convergence 
should occur with fewer modes within that range of ka. 


FMM: Soft Sphere, Aspect Ratio = 1:1, Digits=16 



Figure 7. Original FMM and improved FMM for »/=l, n=\. This is the first plot 
of results achieved with the improved model. With the improved FMM 
calculated at digits=\6, which is the same as the original FMM, the 
identical curves give confidence that VPA was implemented correctly. 

The improved FMM was also run for a prolate spheroid with AR=5:\. The results 
for the spheroid are shown in Figure 9. The improved FMM computed modes m=\, n=\ 
and m=\, n=2 for all four values of precision, but no further modes could be computed 
because of errors. Again, the original FMM results were duplicated by the improved 
FMM for the modal combination m=\, n=\. 

For bothy4i?=l:l andy4i?=5;l, the results of the improved FMM were identical for 
mode m=\, n=\ and also for m=\, n=2 for all four values of digits (i.e., 16, 24, 32, and 
40). This result, which is shown in Figure 10, suggests that SVD is affecting the solution. 
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Figure 8, 


FMM: Soft Sphere, Aspect Ratio = 1:1, Digits=32 



FMM: Soft Sphere, Aspect Ratio = 1:1, Digits=32 



FMM: Soft Sphere, Aspect Ratio = 1:1, Digits=32 



Modes ffi=l, «=1; m=l, n=2; and m=l, n=3 for the original FMM and 
improved FMM with digits=32 for a soft sphere. At digits=32, the 
improved FMM runs at twice the precision of the original FMM. 
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FMM: Soft Prolate Spheroid, Aspect Ratio = 5:1, Digits=32 



FMM: Soft Prolate Spheroid, Aspect Ratio = 5:1, Digits=32 



Figure 9, Modes m=\, n=l and m=l, n=2 for the original FMM and improved 
FMM with digits=32 for a prolate spheroid with ^i?=5:l. These plots 
show distinct differences between the two models for these modes. The 
identical curves in window (a) give confidence that VPA and other 
techniques were implemented correctly. Precision of digits=32 is 
approximately double the precision of the original FMM. 
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FMM: Soft Prolate Spheroid, Aspect Ratio = 5:1, Digits=16 FMM: Soft Prolate Spheroid, Aspect Ratio = 5:1, Digits=24 



FMM: Soft Prolate Spheroid, Aspect Ratio = 5:1, Digits=32 FMM: Soft Prolate Spheroid, Aspect Ratio = 5:1, Digits=40 



Figure 10. Mode »i=l, n=2 of the improved FMM for a prolate spheroid with 
AR=S'.\ at four different levels of precision— digits=\6y 24, 32 and 40. 
The plots are identical, suggesting that SVD or another factor is 
affecting the results. 

C. ANALYSIS AND DISCUSSION 

1. Gains in Converged Range of ka 

For the purposes of analysis, the eonverged ka region of the FMM is defined as 

the value of ka under whieh two eonsecutive modes overlap. The progression of the 

eonvergenee of the original FMM is shown in Figure 11, with the three initial modes 

plotted on top of eaeh other. With the addition of eaeh mode, the eurves overlap eaeh 

other at higher frequeney values. Figure 11 illustrates that the original FMM is 

eonverged for modes m=\, n=2 and m=\, n=3 out to ka=\.2. The progression of 

eonvergenee for the three initial modes is also plotted for the improved FMM in Figure 

12. Notiee that the two initial modes (i.e., m=\, n=\ and m=\, n=2) do not overlap in the 

improved FMM solution, although they may be eonverged somewhere in the ka range 
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that is below the plotted range. However, the next two consecutive modes (i.e., m=\, n=2 
and m=\, n=3) are converged in the illustrated frequency range. 

Figure 11 and Figure 12 can be compared to show that modes m=\, n=2 and m=\, 
n=3 are converged to ka=\.2. This region is expanded in Figure 13 to show a close-in 
view of the point of convergence in both the original and improved FMM. For modes 
m=\, n=2 and m=\, n=3, Figure 13 confirms that the solution converges out to ka=\.2. 


Original FMM: Soft Sphere, Aspect Ratio = 1:1 



Figure 11. Original FMM progression of convergence with addition of higher 
modes. The point of convergence is extended to higher values of ka as 
higher modal combinations are added to the solution. The first null, 
caused by computed scattering interference patterns, is shown by the 
arrow at ka=A for mode m=\, n=3. 
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Improved FMM: Soft Sphere, Aspect Ratio = 1:1, Digits=32 



Figure 12. Improved FMM progression of convergence with addition of higher 
modes. The point of convergence is extended to higher values of ka as 
higher modal comhinations are added to the solution. The improved 
FMM exhibits as much stability as the original FMM out to A:a=10. The 
first null is reached in mode m=l, n=\ at ka=4.S. The null is shifted to 
the right to ka=6.5 for m=l, n=3, which is 2.5 units higher into the 
frequency range than the original FMM predicted for the same mode. 
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Convergence of the FMM: Soft Sphere, Aspect Ratio = 1:1 
Original FMM and Improved FMM at Digits=32 
for Modes m=1, n=2 and m=1, n=3 



Figure 13. Convergence of the original FMM and improved FMM for modes m=l, 
n=2 and m=l, n=3. The dashed lines indicate the two consecutive modes 
of the original FMM while the solid lines indicate the two consecutive 
modes of the improved FMM. The value of ka at which the two curves 
depart from one another is very close for both models at about ka=1.2. 

While it appears as though the eonverged, or overlapping, portion of the eurves 
has not changed with the model improvements incorporated into the improved FMM, 
there is another feature evident in Figure 11 and Figure 12. The null, which results from 
destructive interference patterns in the various functions of the model calculations, is 
shown at ka=4 in the original FMM for mode m=l, n=3 in Figure 11. Figure 12 shows 
that this null is shifted to a higher frequency of ka=6.5 for the same mode in the improved 
FMM. This increase of 2.5 in the range of ka marks a significant performance gain for 
the improved FMM. This suggests that the overall solution, while not overlapping 
beyond the point at which the original FMM does for these low modes, is a better 
characterization of the scattering phenomena and that it is closer to the true converged 
solution. This concept is depicted more intricately in Figure 14, which shows a 
comparison of the original FMM and improved FMM at mode m=\, n=3 with the original 
FMM converged solution at m=\5, n=30. Considering the reduced target strength (RTS) 
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over the entire range of dimensionless frequeney {ka), rather than just the portion where 
consecutive modes overlap, the improved FMM appears to be a more stable 
approximation to ka=6.5. 


Comparisons for a Soft Sphere, Aspect Ratio = 1:1 
Originai FMM and improved FMM with Digits=32 



Figure 14. Comparison of mode m=\, n=3 of the original FMM and of the 
improved FMM at digits=32 to the converged solution of the original 
FMM. The first null space is reached at ka=3.8 for the original model, 
while the improved FMM shifts the first null to ka=6.5. 

2. Accuracy Gains 

An important finding is that the improved model run at digits=\6 yields a 
different converged curve than the original model does for modes m=l, n=2 and higher, 
although the original model also incorporates double precision, or 16 digits into its 
answers. Mode m=\, n=\ of the improved FMM exactly matches the same mode of the 
original FMM. However, there is a difference in the converged portions of the two 
models at higher modes, shown by an arrow in Figure 12, that illustrates the variation in 
model calculations between modes m=\, n=\ and m=\, n=2 for the improved FMM, 
which was not apparent in the original FMM curves displayed in Figure 11. The 
difference exists because the second term in the improved FMM summation is different 
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than that incorporated into the original FMM. The explanation is that the improved 
model performs caleulations with symbolie variables, which are not rounded as often as 
the double preeision floating-point values throughout the course of model computations. 
The values that are set as initial eonditions are exact in the improved FMM, whereas they 
are rounded approximations in the original FMM. Therefore, the improved FMM is more 
accurate even with digits=\6 —^the typical number of digits in the output of double 
precision calculations in the original FMM. 

The difference noted above can be seen more closely in the converged curves of 
the original and improved FMM in Figure 13, whieh are separated by 3.4 dB along the 
RTS axis. A 3.4 dB gain in accuracy for that range of ka values eould be another 
signifieant performance gain for the improved FMM. However, this accuracy difference 
cannot be confirmed until a larger number of modes can be calculated with the improved 
FMM, after the MATLAB errors are surpassed. 

3. Singular Value Decomposition 

The VPA caused errors at lower modal eombinations in the series solution than 
expected, and this prevented significant investigation of the effects of SVD on the 
performanee of the improved FMM. The SVD algorithm is more effective at the higher 
modes, where the numerieal eontribution to the solution is mueh smaller and potential 
error contributions are mueh higher. Since the improved FMM was significantly limited 
in the number of modes that it could produce, SVD could not be investigated any further 
than the useful range of the improved FMM. 

When working with increased precision and symbolic numbers, the threshold 
incorporated into the SVD teehnique must be sufficiently small to rid the solution of only 
those values that are beyond the eomputer’s capability to represent them. In this 
research, symbolic numbers created with vpa and sym were eonverted from symbolic 
form back to double preeision after being carried as far into the solution as possible. 
Since the symbolie variables were carried all the way through the Bessel functions and 
associated Legendre functions, the only rounding incorporated into the variables before 
they reached the SVD portion of the model was through the integration proeess. With 
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this being the case, there was little instability associated with roundoff error incorporated 
into the model up until the SVD commands. 

Although a larger or smaller number of values may be kept in the solution by the 
SVD in the improved model, the values in the matrices are more accurate. Fewer 
rounding operations are incorporated into the symbolic variables of the improved model, 
so fewer roundoff errors are propagated into the matrices. This means that more of the 
right values are being kept in the solution by the SVD numerical filtering process while 
fewer erroneous values are being allowed to stay and affect the results of the scattering 
simulation. 

The SVD algorithm is the most likely cause for the duplication of results for all 
four values of digits, as shown in Figure 10. The SVD may be limiting the number of 
terms that are included in the solution, even for digits=\6. Thus, increasing the precision 
with a higher value of digits accomplishes nothing, because the gains afforded by 
increased precision are negated by the filtering effect of the SVD algorithm. 

4. Shortfalls of Extended Precision and Numerical Analysis Techniques 

Throughout the course of the research, it was determined that the VPA numbers 
could not be carried through the entire model from start to finish, because of inherent 
errors within MATLAB. As mentioned previously, several functions within MATLAB 
have not been adapted from older source code to handle the capabilities of the Symbolic 
Math Toolbox. Although the VPA values could not be carried through the entire model, 
they were carried through the calculation of the Bessel functions and through the 
associated Legendre functions. If the inputs to these functions from previous sections of 
the model code were created in double precision (i.e., the precision of ordinary MATLAB 
numbers), then these inputs would be stored as rounded approximations. Thus, the model 
would have initiated with erroneous values, and the errors would have only been 
amplified throughout further calculations. However, the VPA values with extended 
precision were entered as initial conditions, and the higher precision was carried through 
the furthest point allowed by the MATLAB software before converting the values back to 
double precision. The output of the Bessel function, Hankel function, and associated 
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Legendre funetion ealeulations had to be eonverted to double preeision in order to earry 
out the max command—a command that could not be circumvented. 

Because of all of the errors encountered with VPA, not enough modal 
combinations could be calculated to complete a thorough evaluation of the new 
performance envelope. Even at digits=\6, impassable errors prevented computations 
beyond the combination m=\ and n=3 iox AR=\:\, and results were stopped at m=\ and 
n=2 for AR=5\\. However, the extended precision results that were collected suggest that 
the performance envelope would be positively affected by the model improvements. The 
solutions for the second modal combination (i.e., m=\ and n=2) for both AR=\:\ and 
AR=5:\ showed that the smooth portions of the RTS vs. ka curves extended to higher 
frequencies for the improved model than they did for the original model. This illustrates 
that destructive interference has less of an effect on the curves of the improved EMM, 
and fewer modes may be required to achieve a converged solution. 

Another limitation of the improved EMM is computational expense. The 
computational time required by the model was higher than anticipated, and these results 
are discussed in the following section. This information is presented separately from the 
aforementioned shortfalls, because computational expense is currently the most 
significant limitation of the improved EMM. 

D. COMPUTATIONAL EXPENSE VERSUS PRECISION ANALYSIS 

The EMM requires a great amount of computational expense, especially with 
increased precision. In order to evaluate the processing time requirements added by 
increased precision, the model was run several times at different levels of increased 
precision with all other input variables held constant. 

One glaring item of interest in the RTS vs. ka plots displayed in the previous 
sections is that the improved EMM was only run for the range 0.5<A:a<10 with a ka 
increment of 0.5, which covers only a portion of the RTS vs. ka window in each plot (i.e., 
the improved model curves appear truncated). Eewer values of ka could be processed 
due to the computational time required by the higher values of precision and the symbolic 
mathematics. The original model was run for the range 0.01<A:a<10 with a ka increment 

of 0.01 providing 50 times the number of data points and while requiring only a small 
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fraction of the computational time that the improved FMM used. The plots generated 
with these ka ranges were used for visual comparison only, and a separate set of model 
runs was completed for the computational expense versus precision analysis. 

To compare model run times for the original and improved FMM, both models 
were run for the range 0.5<A:a<10 with a ka increment of 0.5. This meant that the same 
number of data points would be computed by each model. Figure 15 shows the run times 
for the original FMM, while Figure 16 shows the comparable run times for the improved 
FMM. Note the lengthy runtimes associated with the improved model that are labeled 
along the y-axis in Figure 16 in comparison to the short runtimes of the original FMM 
that are shown in Figure 15. The modes that were calculated are annotated on the figures. 
Mode m=\, n=3 is displayed only for digits=32 in Figure 16 because this modal 
combination encountered errors at the other three precision settings. The times shown in 
Figure 15 and Figure 16 are the times that were required to calculate all modes up unto 
and including the indicated mode (i.e., the time shown for m=\, n=2 is the time that was 
required to calculate both m=\, n=\ and m=\, n=2). 

The difference in computational expense between the two models was extremely 
large. This was mostly due to the differences between the computational expense of 
floating-point arithmetic and symbolic mathematics. To put this in perspective, mode 
m=\, n=\ for the original model took only 0.12% of the computational time that the 
improved model required for the same calculation with digits=\6. Reciprocally, the 
improved model took 81,600% of the computational time of the original model at the 
same mode and precision. 

A clear relationship between precision and runtime could not be established for 
the original versus improved models because of the differences between floating-point 
and symbolic mathematics. However, it is obvious from the plots that the symbolic 
mathematics format of the improved FMM with higher precision requires a much greater 
amount of computational time than standard IEEE double precision format. The 
relationship between levels of precision for symbolic math calculations alone is linear. 
This is indicated in Eigure 16 by the diagonal line along the bars corresponding to 

calculations of mode m=\, n=2 at the different values of digits in the improved EMM. 
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Figure 15. Original FMM computational expense in minutes. The original FMM 
was run solely at double precision, which corresponds to about 16 
decimal digits. Note the relatively short amounts of time required to 
run the modes, which are annotated along the j-axis. 
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Figure 16. Improved FMM computational expense in minutes versus VPA digits. 

The improved FMM was run at four different levels of precision, which 
are shown along the jc-axis here. Note the extremely long computation 
times. 
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V. SUMMARY AND CONCLUSIONS 


A. NOTED IMPROVEMENTS IN THE CONVERGED SOLUTION 

Variable precision arithmetic afforded a more accurate numerical solution for the 
improved FMM. A 3.4 dB difference was noted between the converged curves of the 
original and improved FMM at low frequencies. The converged solution denoted by the 
positions of the smooth portions of the curves and null spaces in the RTS vs. ka plots 
reaches 2.5 units higher into the frequency range for mode m=\, n=3 of the improved 
FMM because of the executed numerical techniques. Limitations on the improved FMM 
included unavoidable errors in MATLAB, fewer working modal combinations, and 
computational expense. While limitations on the number of modes that could be 
calculated by the improved FMM prevented a complete analysis of the new performance 
envelope, the expectations of improvements in model performance held true. Another 
bound has been broken for the FMM, and the improved model offers better scattering 
predictions. 

B. FEASIBILITY FOR IMPLEMENTATION 

The FMM requires many improvements before it would be ready for 
implementation into an operational system. However, the performance envelope 
discussed by Reeder and Stanton (2004) proves that the model is useful for a large range 
of modes before increased precision is even considered. 

The FMM is intended for use in predicting acoustic scattering by complex 
scatterer shapes. In light of the computational expense of the improved FMM when 
describing simple shapes like spheres and prolate spheroids, the computational expense in 
the case of scatterers with complex or irregular boundaries is currently prohibitive to 
operational use. This research shows that the Fourier matching method is still far too 
computationally intensive for current computers to handle as a routine calculation. With 
the use of an active sonar system, operators need results within seconds to a few short 
minutes. The FMM developed by Reeder and Stanton (2004) can take up to several 
hours to run various modes in a single scattering solution, making it impractical for use in 
today’s operational arena. Incorporating increased precision into the formulation and 
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model source code extends the computational time even further. Higher precision 
computing only exacerbates the problem of expense, and significant advances in both 
computer processor speed and storage capacity must be made before higher precision can 
be utilized in operational systems. However, it is important to maintain the perspective 
that the FMM is currently a fundamental, physics-based, research level model that is still 
in developmental stages. The FMM pushes the envelope of technology by demanding 
faster and more robust computer systems before it can be fully implemented. 

C. RECOMMENDATIONS FOR FUTURE WORK 

The first four recommendations in this section are provisional solutions that will 
allow more expedient research of the capabilities of the FMM. The fifth and final 
recommendation here will bring a major change to the mathematical formulation of the 
FMM, which may offer significant technical improvements to the future FMM. 

1. Implementation of Symbolic Mathematics 

The results of symbolic operations are exact, since they use rational arithmetic, so 
this may be an area to consider further research for the FMM. While the sym command 
introduces even more errors than vpa in the current version of the MATLAB software, 
several of these errors will no longer be encountered once Mathworks resolves 
compatibility issues with the Symbolic Math Toolbox. The sym command works 
similarly to the vpa command, so implementation would be a relatively straightforward 
extension of this research. Implementation of the sym command would be beneficial for 
the development of the FMM once the issues of MATLAB errors and computational 
expense are mitigated. 

2. Conversion to Fortran 

Converting the FMM to run in the Fortran programming language could permit 
the use of supercomputers for FMM development. Fortran could enable the employment 
of machines with higher values of machine precision, including computers that are 
currently capable of running with 128-bit word storage. The results of this research show 
that increased precision makes a difference in FMM scattering predictions. The full 
capabilities of the FMM may be realized if it can be run at higher modal combinations 
with increased precision without running into errors associated with software. 
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3. Investigation of SVD Thresholding 

Future versions of the MATLAB software may allow the errors experieneed in 
this researeh to be avoided. This would allow the eomputations of higher modal 
eombinations with the improved FMM. At the higher values of precision afforded by 
VPA, the threshold for the SVD algorithm should be set to reflect the minimal amount of 
roundoff error according to the appropriate value of eps. This latest version of the FMM 
incorporates VPA, which allows various settings of precision. Thus, the range of 
appropriate thresholds for SVD at various levels of computer precision should be 
investigated to a greater extent. The extent to which SVD is affecting the solution is still 
unclear. This is the next step in the progression of FMM development, if the improved 
FMM from this research is intended for future use. 

4. Database of FMM Results 

An improvised solution to the lengthy runtimes of the FMM that would support 
implementation may be to tabulate results in a database. A database of FMM results may 
provide operators with the answers they need in near real-time, rather than having to 
calculate the lengthy FMM solution from start to finish. If a compilation of results can be 
generated for use as a look-up table or computer-accessed database, then application of 
those results to operational problems could possibly be accomplished within a more 
reasonable amount of time. 

5. Incorporation of Spheroidal Wave Functions into the FMM 

In 1957, Carson Flammer published a text called Spheroidal Wave Functions 
(Flammer, 1957). This text presents the foundations for the next major revision in the 
mathematical formulation of the FMM. The purpose of this future improvement is to 
change the mathematical formulation of the model to match the general scatterer shape 
more closely. This may yield computational advantages, because the exact solution given 
by the spheroidal wave functions is closer to the actual scattering phenomena associated 
with a prolate spheroid than is the approximation from the spherical wave functions. 
Hence, fewer modal combinations may be required to yield the same result achieved with 
the spherical wave functions of the current FMM formulation. The resulting model 
should replicate scattering phenomena for prolate spheroids or elongated irregular bodies 
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more accurately and possibly with less computational expense. Figure 17 shows a 
comparison of the Reeder and Stanton (2004) FMM with the exact prolate spheroidal 
solution. 


FMM & Exact PS (10:1) 
Rigid Prolate Spheroid 




0.5 


ka 

Figure 17. FMM solutions of Reeder and Stanton (2004) compared to exact prolate 
spheroidal solutions. Many lines coincide because of the similarities in 
solutions, obscuring the dotted line exact solutions under the solid 
FMM curves. (From Reeder and Stanton, 2004) 

Application of the resulting FMM with incorporated spheroidal wave functions is 
envisioned for advanced torpedo sonar systems. The new formulation is expected to be 
more effective in reducing clutter in the sonar picture of a torpedo using active 
transmissions. Inclusion of spheroidal wave functions would also make the FMM more 
applicable to detection of elongated scattering bodies, such as submarines, gliders, and 
autonomous underwater vehicles (AUVs) (Reeder and Stanton, 2005). 
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D. CONCLUDING REMARKS 

The FMM formulation for scattering by irregular, finite-length bodies of 
revolution developed by Reeder and Stanton (2004) has been adapted for extended 
precision through the use of the variable-precision arithmetic (VPA) instrument of the 
MATLAB Symbolic Mathematics Toolbox. The mathematical formulation of the FMM 
has not changed, but the computer processes by which the solution is calculated have 
been optimized for increased accuracy and the reduction of roundoff error. The improved 
FMM results confirm a more accurate converged solution at lower modal combinations, 
indicating that the performance envelope of the FMM has been improved. Extensions of 
this research beyond current software limitations may show marked improvements in the 
model’s ability to depict real-world acoustic scattering events. 

While implementation of the FMM into operational systems is currently not 
gainful due to operational time constraints and other limitations, this research shows 
improvement in the ability of the FMM to represent scattering phenomena. With 
continued performance gains and functional progress, the FMM will be an integral part of 
a superior active sonar system. 
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